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Abstract 

A constant-time solution of the continuous Global Optimization Problem (GOP) is obtained by 
using an ensemble algorithm. We show that under certain assumptions, the solution can be guar- 
anteed by mapping the GOP onto a discrete unsorted search problem, whereupon Briischweiler's 
ensemble search algorithm is applied. For adequate sensitivities of the measurement technique, the 
query complexity of the ensemble search algorithm depends linearly on the size of the function's 
domain. Advantages and limitations of an eventual NMR implementation are discussed. 
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I. THE GLOBAL OPTIMIZATION PROBLEM 



Optimization problems are ubiquitous and extremely consequential. Their theoretical and 
practical interest has continued to grow from the first recorded instance of Queen Dido's 
problem [l|_to_present day forays into complexity theory and large scale logistics applications 
(see Refs. j^l, S, U, 0, |(| and references therein). 

The formulation of the Global Optimization Problem (GOP) is deceptively simple: find 
the absolute minimum of a given function - called the objective function - over the allowed 
range of its variables. Sometimes, the function whose global minimum is to be found is not 
specified in analytic form and must be evaluated by a computer program. When we can 
only access the output of the computation, the program acts as a black-box tool, called an 
oracle. 

Of course, the brute force approach of evaluating the function on its whole domain is 
either impossible if the variables are continuous, or prohibitively expensive if the variables are 
discrete, but take values in large domains and in high dimensional spaces. Since, in general, 
each oracle invocation (function evaluation) involves an expensive computational process, 
the number of function evaluations needs to be kept to a minimum. It is not surprising 
that, together with accuracy, this number probably provides the paramount criterion in 
comparing the efficiency of competing optimization algorithms. 

The primary difficulty in solving GOPs stems from the fact that the condition for deter- 
mining minima, namely annulment of the gradient of the objective function, is only necessary 
(the function may have another type of critical point) and local i.e. it does not distinguish 
between local and global minima. The generic strategy to find the global minimum involves 
two main operations, namely: (i) descent to a local minimum and (ii) search for a new de- 
scent region. However, this generic strategy is powerless for certain problems. The following 
one-dimensional example (the golf course problem) illustrates the difficulties. Define the 
function / : [0, 1] by 

f for < x < a - e/2 and a + e/2 < x < 1 

M = { ~ ~ ' (i) 

I -1 for a -e/2 < x < a + e/2 

where < e C 1, and a is a point in the interval (e/2, 1 — e/2), but is otherwise unknown. 
From this definition, it is clear that to find the absolute minimum of /, one should evaluate 
it within the e well-interval around the unknown point a. If this function is defined like 



an oracle, i.e. if one does not know its definition above and the value of the number a, 
the probability P e of choosing the variable x within the well-like interval where f(x) = — 1, 
is e. In the n-dimensional version of this problem, the probability P t becomes e n , so the 
complexity of the problem grows exponentially with n (the dimensionality curse). Moreover, 
knowledge about the (partial) derivative(s) of / would not help, since they are all zero 
whenever defined. In the light of this example, it seems that without additional information 
about the structure of the function there is no reasonable expectation of deciding upon 
an efficient optimization strategy and one is left with either strategies that have limited 
applicability or the exhaustive search option. 

Thus, new approaches are needed, which take advantage of additional information to 
reduce the complexity of the problem to a manageable level. We emphasize that this infor- 
mation is actually available in some classes of GOP - for instance it happens to be available 
for the golf course problem above - but cannot be taken advantage of within conventional 
optimization algorithms. 

Here we present an approach that uses this information to map the continuous GOP 
into a discrete search problem (DSP). Once this mapping is completed, we apply a variant 
of Briischweiler's search algorithm, as described in the following section, to obtain a value 
within the basin of attraction of the global minimum in a number of function evaluations 
that scales logarithmically with respect to the cardinality of the resulting discrete set. While 
the parallelism of ensemble computing indeed provides an "exponential speed-up", it also 
demands an exponentially large number of processors. Thus, we use the phrase "constant- 
time" , to account for the relationship between speed and resources. 

II. MAPPING THE CONTINUOUS GOP TO A DISCRETE SEARCH PROBLEM 

Consider a real function of d variables, /(x), x = (x±, X2, ■ ■ ■ , Xd) that has a certain 
degree of smoothness. The precise definition is not critical for the argument below. Without 
restricting generality, we can assume that / is defined on [0, l] d — > [0, 1]. Assume now that: 
(i) there is a unique global minimum and its value is zero; (ii) there are no local minima 
whose value is infinitesimally close to zero; in other words, the values of the other minima 
are larger than a constant 5 > 0, and (iii) the size of the basin of attraction for the global 
minimum measured at height 5 is known to be r^. We note that these assumptions can be 
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actually verified for the problem at hand, if additional information is available. Also, these 
assumptions may in general be relaxed to some extent, as outlined in the first part of the 
Discussion. 

The implementation idea is the following. Instead of /(x), consider the function <?(x) := 
(/(x)) 1 /" 1 . For sufficiently large m, this function will take values very close to one, except 
in the vicinity of the global minimum, which will maintain its original value, namely zero. 
Of course, other transformations may be used to achieve essentially the same result. Now 
divide the hypercube [0, l] d in small (^-dimensional hypercubes with sides Ax = 1/M, where 
M is a positive integer. At the midpoint of each of these hypercubes, define the function 
h{x) := 1 - INT[g(x) + 1/2] (INT denotes the integer part). 

As a result of the transformations above, we reduced the GOP to a DSP for h : 
{1,2, . . . , N} — ► {0,1}, which is an integer-valued function defined on a discrete set of 
N := 2 n+1 points. The function h(i) is known to be zero for all inputs i, except for one 
special input % = q, where h(q) = 1. Here N = M d and the domain on which the function h 
is zero depends on m and on the size, t$, of the basin of attraction of the global minimum 
at height 5. Thus the problem becomes to find the value of the special input q efficiently 
Recently, we proposed an adaptation of Grover's quantum search al gori thm to solve the 
continuous GOP 0, which requires 0(y/~N) function evaluations P. B . llo|. 

In this paper, we propose to further improve the efficiency, by using an ensemble search 
algorithm, which is described in detail in the following section. The algorithm is adapted 
from Bruschweiler's ensemble search algorithm 11. l^. using a parallel technique to achieve 



exponential speedups for problems up to a certain size, dictated by the sensitivity of the 
measurement process. As mentioned before, efficiency is understood in relation to the query 
complexity of the algorithm. Indeed, when N is large and the function evaluation (i.e. the 
oracle) is costly in terms of computational complexity, reducing the number of function 
evaluations is critical. 

Application of the ensemble search algorithm to the function h results in a point that 
returns the value zero. It is easy to see that, by construction, this point belongs to the 
basin of attraction of the global minimum. We return then to the original function / and 
apply the descent technique of choice that will lead to the global minimum. If the basin of 
attraction of the global minimum is narrow, the gradients of the function / may reach very 
large values which may cause overshots. Once that phase of the algorithm is reached, one can 
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apply a scaling (dilation) transformation that maintains the descent mode but moderates 
the gradients. On the other hand, as one approaches the global minimum, the gradients 
become very sma . 
may be required 

the algorithm above will find it without additional complications. 



1 and certain acceleration techniques based on non-Lipschitzian dynamics 
lij^l- If the global minimum is attained at the boundary of the domain, 



III. BRUSCHWEILER'S ENSEMBLE SEARCH ALGORITHM 



Following Grover and Shor's breakthrough algorithms for quantum computing [8|, [13|, 3| > 
an alternative paradigm for computing has been suggested by Madi, Briischweiler, and 
Ernst, which operates on ensembles i.e., mixed states of identical spin sytems, using a spin 
Liouville space formalism for density operators The Liouville space for a system of n 
spins contains the n x n density matrices representing all the possible quantum states of 
the spins. Ensemble algorithms are not quantum algorithms, insofar as they do not involve 
entanglement of quantum states, so they only use a small subset of the Liouville space. 
Throughout this paper, "mixed states" describe a statistical ensemble of spins, rather than 
individual spins, so that each element of the ensemble performs part of the computation, in 
the same way as each processor in a classical parallel computer. 

This new paradigm exploits the parallelism offered by simultaneously acting on many 
different input states in an ensemble of spins. In contrast, quantum computing with pure 
states relies on the parallelism of entangled states, to perform operations in a Hilbert space 
which is the tensor product of multiple qubits. 

The advantage of ensemble algorithms is that they may be exponentially faster than the 
equivalent quantum algorithms, for adequate measurement sensitivities. Two realizations of 
ensemble algorithms are provided b y B ruschweiler's search algorithm and the summing 
algorithm proposed by the authors [lfij . which illustrate the trade-off between memory and 
speed capabilities. 

Bruschweiler's algorithm for searching an unsorted database, uses the ensemble computing 
paradigm in the context of NMR technology. The strategy employs binary partition of the 
N elements in the database, to find the desired element after 0(log 2 N) oracle queries. We 
note that the lower bound derived by Nayak and Wu for the query complexity of quantum 
search algorithms, does not apply here, since the quantum lower bound used in their 
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derivation assumes that the initial state is a pure state. 

Ensemble computing requires an exponentially larger set of physical resources to encode 
the same number of distinct input states compared to quantum computing with pure states. 
On the other hand, ensemble computing holds the important advantage that it is insensitive 
to the decoherence time of the spins, which is an outstanding limiting factor for quantum 
computations involving entangled states. 

A "divide and conquer" scheme is used to test whether the special input q belongs to 
exponentially finer and finer partitions of the set of input values. The algorithm is envisaged 
for a physically realizable system consisting of binary-valued spins, therefore the number of 
spins needed in the input register varies from n to 1, to accommodate the number of input 
values in partitions ranging in size from N/2 values down to 2 values, respectively. It is also 
assumed that a 1-spin output register is available, to encode the value of h e {0, 1}. 

The first step of the algorithm is to divide the input values into two equal partitions, 
X\ = {1, . . . , N/2} and Xi — {N/2 + 1,...,N}, and test whether the special input q 
belongs to Xi or Xi ■ The numerical subscript, k — 1, indicates the size of the partition as 
a fraction (l/2 fc ) of the entire set of input values, x = {!>••• > an d the superscript ± 
differentiates between the lower and upper partitions. 

The n-spin input register is initialized as an ensemble, namely an equally-weighted mixed 
state, 

Pin = JjJ^ E K)n(i|n, (2) 

which accounts for half of the N = 2 n+1 possible input configurations. The density operator 
can be represented in Liouville space by a density matrix that has non-zero elements only 
on its diagonal. The off-diagonal elements are all zero, indicating the absence of quantum 
coherence between any of the states \i) n . The ket states \i) n correspond to the eigenstates 
of the Zeeman Hamiltonian created by a strong external magnetic field , and can also be 
written in terms of individual spins, 

\i) n = | a»i) <8> \a i2 ) ® ■ ■ ■ ® \a in ) (3) 

where (an, a^, . . . , a in ) e {0, 1} are the digits of the number (i — 1) in binary format, and 
g) denotes the tensor product between the states. The bra states, (i\ n , are the adjoint of 
the ket states. For further convenience, the state |0) denotes a spin "up" and the state |1) 
denotes a spin "down". 
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Each of the N/2 input configurations in the partition \i is realized by a sub-ensemble 
which contains many identical molecules with n spins. The number of molecules in each 
sub-ensemble may vary, but the higher it is the more robust the algorithm is going to be 
against random errors. Ideally, all of the molecules in each sub-ensemble are initially set to 
one of the input configurations \i) n (i\ n , and the entire sub-ensemble responds to external 
controls as one "mega-spin". 

At room temperature, the thermal equilibrium state of an ra-spin ensemble in an NMR 
experiment is equal to the sum of p\^ and a traceless density matrix, with zero off-diagonal 
terms Q| • The diagonal elements of the traceless density matrix are given by permutations 
of Boltzmann factors, which have small magnitudes with respect to 2/N, therefore the 
thermal equilibrium state closely approximates the desired initial mixed state p^. 

The original approach proposed in Briischweiler's search algorithm requires that the out- 
put register is initially set to zero. This requirement cannot be satisfied in an NMR imple- 
mentation, as it is impossible to set all the spins in the output register to |0). However, 
Bruschweiler's search algorithm is still applicable, by taking the initial state to be the natu- 
ral state available in NMR: a thermal mixture of |0) and |1) spins. Thus the entire ensemble 
(input and output registers) is described by p^ ® Pout, 

f$ ® P»* = jjfi Y, N)n(i|n®^{|0)(0| + |l)(l|}. (4) 

The effect of this modification is that measurements of the output register need to distin- 
guish deviations from a thermal mixture, rather than from a uniform ensemble of |0) spins, 
following application of the search algorithm. 

The function h, analogous to the oracle in Grover's search algorithm, is evaluated by 
applying a reversible unitary transformation Uh defined as: 

U h \i) n ®\0) -> |i) n <g> \h(i)) (5) 
U h \i) n ®\l) -> |i)„®|l). (6) 

The transformation has no effect on the eigenstates \i) n , but partitions the output register 
into two sub-ensembles: 1/2 — 1/iVof the output spins have the value |0), and 1/2 + 1/N 
of the output spins have the value |1). The action of Uh on the initial mixed state can be 
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written as 

u h (p^®Pout)ut = j±- J] |<>n(<U®i{IMOXMOI + |i><i|}- (7) 

The transformation [7/j can be described jl^] as a permutation of all states spanned by 
the input and output registers, which is applied to the system, to evaluate the function 
h simultaneously on all of the input values % given by the initial state in Eq.(j3J). This is 
equivalent to evaluating the function h concurrently on N/2 processors in a classical parallel 
computer. In a physical implementation, Uh is the product of a sequence of fundamental 
unitary transformations for each of the spins in the output register. 

The effect of the transformation Uh is to create a small imbalance in the initial thermal 
state of the output register. If we can detect this imbalance reliably, then the special input 
value we are searching for is contained in the partition Xi ■ 

In the physical implementation of the algorithm, the output register generates an output 
signal proportional to the number of spins in state |1). Hence the output signal will increase, 
if xl does contain q, or it will remain the same, if Xi does not contain q. The latter case 
implies that xl must contain q. This ensemble measurement is an analog process, which 
performs a single concurrent evaluation, and serves as a definite test of whether q belongs 
to Xi or Xi- Since only one out of N/2 input values may change the output register, 
under the transformation Uh, the corresponding increase in the output signal is rather small 
(proportional to 2/N). 

Based upon this measurement, the partition containing q is selected and then subdivided 
into two equal partitions, xi an d xt- The equally- weighted mixed state p 4 -™ ^ corresponding 
to the input values in x~2 is prepared in an (n — l)-spin input register, 

pt~ l) = Jfjl Yl N)n<i|n, (8) 

and the same approach outlined above is followed, to test if the special input q belongs to 
X2 or xt- We note that, from an experimental point of view, it could be more convenient 
to prepare the same equally-weighted mixed state in the input register, regardless of which 
partition actually contained q, and then to use this information to determine an equivalent 
transformation that realizes the oracle for the appropriate partition. 

Each test to determine the partition that contains q, reduces the size of the set of input 
values that contains q by a factor of 2, therefore the special input value will be determined 
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by a series of n tests, i.e. the query complexity of the algorithm is 0(log 2 N). The final test 
involves two equal partitions, Xn an d Xn> each containing a single input value, one of which 
is q. 

For the first partition test, it is necessary that the measurement be able to distinguish 
changes in the output signal with a precision equal to or better than 1/JV. This requirement 
may be gradually relaxed for subsequent tests. Indeed, at the k-th test, the precision of the 
measurement only needs to be equal to or better than l/2 fc , and by the final test, which 
determines whether q belongs to Xn 01 Xn> the precision has to be equal to or better than 
1/2. Of course, setting the precision to the minimum required level, l/2 fc , for the k-th 
partition test, may provide some practical advantages in an experimental implementation, 
though it does not lower the query complexity. 

An additional limitation on the measurement sensitivity for the output register is the 
inevitable noise present in measurements, especially due to thermal fluctuations in an NMR 
implementation. As the number of sample points iV = 2 ra+1 increases, the difference be- 
tween distinct normalized output signals in the first partition test will eventually become 
comparable to the noise level. In order to reliably measure the maximum normalized output 
signal l/2 ra in the first partition test, we require an error level Ai < l/2 n+1 . If the error 
level for a single experimental trial is Ai > l/2 fe {k < n), we cannot reliably measure the 
signal from the output register for any of the first k — 1 partition tests, because significant 
differences, which are smaller than l/2 fc , cannot be detected reliably. 

However the error can be reduced by repeating the algorithm a number of times. The 
error level scales inversely with the square-root of the number (N e ) of experimental trials, 
A]y e oc l/y/Ne, because the variance in the measurement is inversely proportional to N e . 
This implies that to successfully implement the first partition test, we have to repeat it at 
least N e PS 2 2<n+1 ~ fe ) = (iVAx) 2 times. Thus the single-run query complexity of the algorithm 
is 0(log 2 iV), but the overall query complexity becomes 0(iV e log 2 iV), due to the repeated 
trials. When k ~ n, N e ps 0(1), so the query complexity remains 0(\og 2 N). However, in 
the worst case, k <C n, the query complexity increases to 0(N 2 log 2 N). 
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IV. DISCUSSION 



The Optimization Problem. We have shown that a broad class of GOPs can be solved 
in a constant time, with a guarantee, by exploiting available additional information about 
the objective function. Our approach enables us to map the GOP into a discrete search, 
and guarantees to find a solution for the global minimum in a constant time, by efficiently 
eliminating all basins of attraction except for the basin of attraction that contains the global 
minimum. 

The additional information obtained by satisfying the assumptions (i)-(iii) in Section HTl 
is readily available for some notoriously hard GOPs, such as the golf course problem. In 
general, a certain amount of information is usually available in every application, but present 
algorithms either do not or cannot use it. Sometimes, this information may be difficult 
to obtain, but its benefits significantly outweigh its cost. Indeed, for high- dimensional, 
computationally intensive problems, the advantage of reducing the complexity of the problem 
from exponential to polynomial is extremely significant. On the other hand, if no information 
is available at all, then a solution for the most general type of continuous GOP cannot be 
guaranteed in finite time. 

In principle, the conditions (i)-(iii) imposed on the functions / in Section |Il] are satisfied 
can be relaxed to some extent. 

Assumption (i) is satisfied by a large class of important practical problems, namely pa- 
rameter identification, pattern recognition, and other classes of inverse problems. In these 
problems, the absolute minimum, namely zero, is attained for the correct values of the 
parameters, matching of patterns, and fitting of output to input. Assumption (i) can be 
relaxed in the sense that the function may have multiple global minima, all equal to zero. 
Functions with multiple global minima will simply result in search problems with multiple 
"special" elements and can be treated accordingly. 

Assumption (ii) can be replaced with the much weaker and altogether reasonable assump- 
tion that / has a finite number of local minima. This would prevent the value of any local 
minimum to be infinitesimally close to the value of the global minimum. 

Assumption (iii) is the most difficult to relax in practice, since this information is crucial 
in ensuring the high efficiency and guaranteed success of the algorithm. In fact, without 
additional information about the size of the basin of attraction, there is no algorithm that 
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- applied on a general problem - can guarantee to find the global minimum in reasonable 
time. Indeed, if the discretization cell is much smaller than the size r$ of the basin of 
attraction, many input values will return the output zero upon the application of the oracle 
Uh- In order to determine one of these input values with certainty, the ensemble search 
algorithm has to be applied for the same number of steps used to find a single special input 
value. Any of the determined values belong to the basin of attraction and will eventually 
lead to the global minimum. On the other hand, if the discretization cell is larger than 
r$, there is a chance to miss the basin of attraction altogether. Moreover, if the basin of 
attraction is infinitesimally narrow, any algorithm would require an unbounded number of 
evaluations. 



The Search Algorithm. The overall query complexity of the ensemble search algorithm 
used in this paper is always less than the overall query complexity for Grover's search algo- 
rithm implemented using NMR pseudopure states [17| . which is 0(N 2 \/N). However, the 
ensemble search algorithm is less efficient than the (theoretical) implementation of Grover's 
search algorithm using pure states, for database sizes above a certain threshold value, which 



is determined by the error level. Briischweiler 12j estimates that his algorithm is more 
efficient only for databases of size N fulfilling the condition 

Ny/Nlog 2 N < A^ 2 . (9) 

We note that it is appropriate to make the above comparisons, since adding more quantum 
processors in an attempt to parallelize Grover's search method, does not achieve any further 
reduction in complexity. 

In addition, ensemble algorithms do not require quantum coherence, thus they are more 
robust and demand less additional resources than quantum algorithms. For example, in 
an NMR implementation, quantum algorithms need to finish within the T 2 decoherence 
time of the system. This requirement is much more stringent than that for NMR ensemble 
algorithms, which only need to finish within the thermal relaxation time T\ m 50-1000 T 2 

00. 

To date, Briischweiler's search algorithm has been implemented experimentally by two 
separate eroupsH Q, b0th using a 3-qubit NMR ensemble system. While a database 
size up to N = 8 is not very impressive, these experiments have demonstrated that the 
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algorithm is based on sound principles. Both groups have improved on the original form of 
the algorithm, by either using measured spectra [20] to read out the ensemble average of 
the output register, or by completely removing the need for the output register [2l|, and 
using a modified oracle transformation, so that the algorithm becomes irreversible. These 
modifications make the algorithm more robust, and easier to realize, respectively, but they 
do not change the query complexity results we presented in Section ITTI1 Other searches in 
NMR ensemble computing have also been proposed and implemented ^'l. . 

However, we note that a number of important experimental issues remain to be addressed, 
before the algorithm can be scaled up to very large numbers of spins. These issues include 
finding suitable molecules with n > 3 nuclear spins, designing the sequence of pulses that 
implements the oracle transformation Uh efficiently, given that the interaction of remote 
nuclear spins is very weak, and enhancing the measurement sensitivity for large ensemble 
systems. The point of this proposal is proof-of-principle, which may eventually lead to 
something practical once NMR technology develops significantly. 

The SNR for a single measurement is currently 10 4 — 10 7 depending on the type of nuclei 
used (24}], and the strength of the applied magnetic field. Taking the best corresponding 
value for the error level, A x « 1CT 7 , the ensemble search algorithm is more efficient than 
Grover's search algorithm using pure states, for database sizes up to iV 10 8 , which 
corresponds to an input register of n m 27 qubits. This threshold size is already of practical 
interest, and will increase in the future, as the error level decreases with advances in NMR 
technology. The enhancement of measurement sensitivity in NMR implementations is 
presently being investigated by a large number of research groups, which are aiming to 
decrease the error level by some constant factor. In the future , h igh-sensitivity methods 
from other fields such as optical sensing of ultra-weak signals [25] using injection-locked 
lasers may also be adapted to NMR experiments. 



V. CONCLUSION 

In this letter, we have shown that additional information about the objective function can 
be used to reduce the complexity of the GOP to manageable levels, thus guaranteeing that 
the correct solution will be found. Conventional classical algorithms cannot take advantage 
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of this information and quantum algorithms yield only a polynomial reduction in complexity. 
In contrast, ensemble algorithms hold the prospect of solving the continuous GOP in a 
constant time. For a restricted number of function samples, where the database size satisfies 
Eq.p|). an ensemble search algorithm requires only 0(\og 2 N) invocations of the objective 
function /. However, this result is not asymptotic with respect to the input size of the GOP, 
due to the limitation imposed by Eq.flHJ). 
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